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The X-56A aircraft is a remotely-piloted aircraft with flutter modes intentionally 
designed into the flight envelope. The X-56A program must demonstrate flight control while 
suppressing all unstable modes. A previous X-56A model study demonstrated a 
distributed-sensing-based active shape and active flutter suppression controller. The 
controller relies on an estimator which is sensitive to bias. This estimator is improved herein, 
and a real-time robust estimator is derived and demonstrated on 1530 fiber optic sensors. It 
is shown in simulation that the estimator can simultaneously reject 230 worst-case fiber optic 
sensor failures automatically. These sensor failures include locations with high leverage 
(or importance). To reduce the impact of leverage outliers, concentration based on a 
Mahalanobis trim criterion is introduced. A redescending M-estimator with Tukey bisquare 
weights is used to improve location and dispersion estimates within each concentration step 
in the presence of asymmetry (or leverage). A dynamic simulation is used to compare the 
concentrated robust estimator to a state-of-the-art real-time robust multivariate estimator. 
The estimators support a previously-derived mu-optimal shape controller. It is found that 
during the failure scenario, the concentrated modal estimator keeps the system stable. 
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Nomenclature 

maximum desired strain variation on sensors upstream of the fiber break 

bias on k th sensor near the fiber break 

current M-step 

number of M-steps 

current concentration step 

number of concentration steps 

squared Mahalanobis distance 

squared Mahalanobis distance of the argument 

upper bound of D 2 

deformations defined at time r 

reference deformations 

finite residuals of all sensors 

finite residual of the k th sensor 

plant 

tuning constant for weight function of M-estimator 

controller 

sensor station 
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estimated modal displacements 
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modal displacement references 
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set of all sensors 
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set of sensors upstream and near the fiber break 
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reference airframe states 

estimated modal displacement states 

Cartesian coordinate in the y-direction 

Cartesian coordinate in the z-direction 

change in aircraft velocity 

change in angle of attack 

change in pitch angle 

population variance-covariance matrix 

k th measurement normal error distribution 

rigid body pitch angle, deg 

arbitrary increase in distribution of squared Mahalanobis distance as a function of argument 
mean of normal error distribution 
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standard deviation of a normal error distribution 

Ofc 

median absolute deviation (MAD) 
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discrete time step 
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deformation modal matrix, a collection of mode shapes. 
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rigid body bank angle, deg 
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strain matrix defined at SFOS measurement locations 
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List of Acronyms 

active flutter suppression 

ASC 

= 

active shape control 

AW IB 

= 

asymmetric wing first bending 

AW IT 

= 

asymmetric wing first torsion 

CME 

= 

concentrated modal estimator 

CPU 

= 

please define, used in figure 7 

FOS 

= 

fiber optic sensors 

IRLS 

= 

iterative recursive least squares 

LTS 

= 

least trimmed squares 

LWLE 

= 

left-wing leading edge 

LWTE 

— 

left-wing trailing edge 

M 

= 

Maximum Likelihood Estimator 

MAD 

= 

median absolute deviation 

MCS 

= 

Monte Carlo simulation 

MM 

= 

Modified Maximum Likelihood Estimator 

NASA 

= 

National Aeronautics and Space Administration 

OLS 

= 

ordinary least squares 

RWLE 

= 

right-wing leading edge 

RWTE 

= 

right-wing trailing edge 

SFOS 

= 

simulated fiber optic sensors 

SW1B 

= 

symmetric wing first bending 

SW1T 

= 

symmetric wing first torsion 


I. Introduction 

T he primary objective of the X-56A (Lockheed Martin, Bethesda, Maryland) program is to demonstrate active 
flutter suppression (AFS) 1 . The experimental flight controllers must suppress flutter modes that have been 
designed into the flight envelope. The long-term goal of the X-56A program is to support extremely lightweight 
flexible structure designs for fuel-bum-efficient aircraft. Lightweight flexible structures may require active control 
to mitigate unfavorable aero-structural coupling. 2 

Part of the X-56A program includes experimental applications of fiber optic sensors (FOS) with fiber Bragg 
gratings in the control system. The FOS measure high-density strain measurements along the entire wing span, ft has 
been shown that simulated fiber optic sensors (SFOS) enable both AFS and active shape control (ASC) for a 
simulated wing model 3 and the X-56A model. 2 

At a certain flight speed, the X-56A models are subject to flutter in the flight envelope. The control system must 
therefore consistently function to ensure the safety of the vehicle. The FOS is part of a sensor suite supporting the 
structural estimation component for the control system, thus, the control system must be tolerant to failures in 
the FOS. 
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Sensor failures must always be expected and prepared for. As such, safety measures will be taken before 
experimental testing of the FOS past flutter speed. For example, if a break in the fiber occurs, the control system 
must continue to function, otherwise flutter could escalate and the aircraft could be damaged or destroyed. 

Researchers at the National Aeronautics and Space Administration (NASA) Langley Research Center (LaRC) 
demonstrated that a break in a fiber 4 * can produce large biases in all downstream sensor measurements. They showed 
that downstream sensors continue to feed back strain, albeit strongly biased strain data. Without proper precautions, 
the control system will continue to utilize the biased data. 

If the control system responds to strongly biased strain, control-induced instability can result. Worse, the control 
system could contribute to the growth of flutter. Therefore, any estimation system that relies on the FOS system 
must be robust to sensor failures through either passive or active means. 

The X-56A simulation model utilizes modal estimates in a previously-developed AFS and ASC controller, 2 
relying on a modal filter to predict the modal displacement states of the vehicle. The modal displacement states are 
fed directly to the controller. The modal filter calculates modal coordinates at every discrete time step by performing 
an ordinary least squares (OLS) on the measured strain, where the column space is the strain mode matrix. The OLS 
has a breakdown point of 0, meaning that even one biased measurement can bias the OLS estimates or modal 
displacements. 

The failure in the FOS (see Ref. 4) could lead to substantially large gross outliers in the sensor data. The OLS 
modal filter thus must be replaced by a practical robust modal filter. This paper presents this solution by introducing 
a new estimator which meets the requirements for a multivariate on-line robust estimation. This estimator supports a 
practical distributed-sensing-based control system. A brief history is presented to provide context for the robust 
modal filter. 

A. Background 

The X-56A program is a joint effort between Lockheed Martin and the United States Air Force Research 
Laboratory to design and develop high-altitude, subsonic, long-endurance autonomous aircraft. 1 5 The aircraft will be 
delivered to the NASA Armstrong Flight Research Center (AFRC) for further experimental research. The 
finite-element models were delivered to AFRC by Lockheed Martin. The models were used to generate plant models 
with aeroservoelastic interactions using the software tool ZAERO. 6 An artistic depiction of the X-56A in flight is 
presented in Fig. 1 (source http://www.wpafb.af.mil/news/story.asp?id=123291532, cited April 29, 2014). 



Figure 1. The X-56A aircraft depicted in flight 
(source http://www.wpafb.af.mil/news/story.asp7kN123291532, cited May 6, 2014). 
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The aircraft comes with a rigid center body and detachable flexible wings; the aircraft is equipped with ten 
trailing-edge control surfaces. All surfaces are available for AFS, ASC, and flight control, although some 
partitioning of duties may be assigned. The aircraft is expected to be tested at subsonic speeds and at low altitude. 2 

1. Fiber Optic Sensor Placement 

Previous computational work demonstrated that shape control was feasible with SFOS feedback in the 
controller. 2 1 In one study, the SFOS was laid out and simulated on the left wing and the right wing of the X-56A 
model. The sensors represent six fibers. Each fiber contains hundreds of strain measurements. Each measurement is 
spaced one-half-inch from the next. The sensor configuration is presented in Fig. 2. 
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Figure 2. The fiber optic sensor layout on the X-56A model. 

The sensor locations shown in Fig. 2 are used to form the SFOS strain mode matrix, T F0S , as described in Ref. 2. 
The points selected for deformation control are located at the right-wing trailing edge (RWTE), right-wing leading 
edge (RWLE), left-wing trailing edge (LWTE), and left-wing leading edge (LWLE). These points were selected to 
maximize modal information for the first symmetric bending and torsion modes. 1 The virtual deformation controller 
tracks modally transformed references, thereby indirectly tracking deformations at these points. 2 

2. Strain Mode Matrix Development and Use for Shape Control 

To simulate strain measurements (or SFOS) measurements, the strain mode shapes computed at the sensor 
locations are given in Fig. 3. 
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Figure 3. Sensor strain mode shapes: a) symmetric wing first bending; b) asymmetric wing first bending; c) symmetric 
wing first torsion; and d) asymmetric wing first torsion. 

The strain mode shapes are originally derived from an MSC Nastran (MSC Software Corporation, Santa Ana, 
California) modal analysis and an elemental strain conversion algorithm presented in Ref. 2. The relationship of the 
measured strains and the strain mode matrix is used to estimate the strain information for SFOS as shown in Eq. (1): 

SkiXc.yc.Zc.T) = Vfos<7(t) + w (1) 


where w is sampled from a noise distribution, q{x) are modal displacement states. The modal displacement states 
can be estimated by a typical OLS modal filter at discrete time r as shown in Eq. (2): 

q(r) = ^FOS s k {x c , y c , z c , r) (2) 


where f is the Moore -Penrose Generalized Inverse 7 and s k (x c ,y c ,z c ,-t) is the SFOS strain measurements. The 
states can then be tracked in the controller. To determine what modal displacement states must be tracked, a modal 
reference is formed from desired displacements on the wing. To accomplish this, the displacement modal matrix 
from MSC Nastran was corrected to be pure elastic and was then used to convert deformation references d re f(r) at 
l r indices to modal displacement references q re f (t) in the manner shown in Eq. (3) 2 : 


Qref(j-') ^ (/ r , m r )dref(T) 


( 3 ) 


where m r is the index of tracked modal displacements in the flight controller corresponding to the symmetric first 
wing bending (SW1B) and symmetric first wing torsion (SW1T) modal displacements. The controller is designed to 
minimize q re f(r ) — q(j), which in turn minimizes d re f(j) — d(r), if modes dominating the response are included 
in q(r ). This process is also referred to as virtual deformation control. The next section presents an overview of 
some statistical characteristics of the strain mode matrix, which make robust estimation of modal displacements 
difficult. 
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3. Statistical Characteristics of the Sensor Strain Modal Matrix 

As previously stated, the purpose of this paper is to robustly estimate modal displacements q(j) in order to 
accurately track q re /(T) even during sensor failures. The problem is multivariate for the X-56A model, because 
14 modes are modeled in the state space model and in the sensor strain modal matrix. Since most theory-based 
methods for robust multivariate estimation assume that unbiased data assume a nominal multivariate normal 
distribution, the test for normality of the sensor strain modal matrix is required. 

The Q-Q plot is a tool that is used to verify if one distribution is similar to another. Here, it is used to verify if a 
sensor strain matrix distribution matches a multivariate normal distribution. If the distributions are similar, then the 
Q-Q plot will result in a line. If the distributions are dissimilar, then the Q-Q line will exhibit unusual behavior in a 
particular direction. The squared Mahalanobis distance is a scale-invariant distance used for one axis of the Q-Q test 
and is given in Eq. (4): 


D 2 = (X- l i m )Z m - 1 (X- l i m y (4) 

where X is the data matrix, ji m is the coordinate-wise population location, £ m is the population variance-covariance. 
It has been shown that the distribution of squared Mahalanobis distance of multivariate normal data assumes a 
chi-square distribution. 7 8 A plot of squared Mahalanobis distance for the strain modal matrix discussed above against 
the quantiles of a chi-square distribution is shown in Fig. 4. 



Chi-square quantile 

Figure 4. Squared Mahalanobis distance versus chi-square quantile for sensor strain modal matrix 

The plot of squared Mahalanobis distance in Fig. 4 skews to the right and then curves strongly upward. The skew 
of the squared Mahalanobis distance is an indicator that the sensor strain data matrix is not multivariate-normal. 
Mardia’s skew and kurtosis estimates 9 indicate that the distribution is subject to large multivariate skew and 
kurtosis, which is non-normal. This is also true when analyzing the multivariate skew and kurtosis corrected for 
small samples. Therefore, it can be expected that some sensors will be much more important than others, indicating 
the presence of leverage sensors or sensors with high importance. 

Since the measured strain is approximately a linear combination of the sensor strain data matrix, the amplitudes 
of the measured strain will also be predictably asymmetrically distributed. Loading will vary with aerodynamic 
condition, thus, the underlying strain distribution will be also difficult to predict. Without a known nominal 
distribution, the application of most computationally efficient robust outlier detection methods is challenging. The 
sensor strain data matrix provides a rough prediction of the nature of strain variation, and may aid in the detection of 
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outliers, since the strain is not expected to vary strongly away from a linear combination of the mode shapes. The 
next section overviews the challenges that are associated with designing a robust modal estimator. 

B. Challenges for On-line Robust Modal Displacement Estimation 

The challenges are many for the robust modal displacement estimator. The estimator must be robust to small, 
medium, and gross outliers. Any estimator which, from a pool of measurements, can resist up to 50% of the 
measurements known to be outliers is known as a high breakdown estimator. 10 The breakdown point is typically 
defined as the percentage of outliers which an estimator can handle before its estimate is biased. A high breakdown 
point is certainly a goal for a robust modal filter. 

The estimator must be robust to at least a single fiber break, which could occur at any time during operations. 
A failure could potentially introduce hundreds of biased sensor measurements simultaneously, which the estimator 
must reject. Depending on how many fibers are used and where the break occurs, the effect of the loss of a fiber on 
the percent of sensors failed will change. 

The controller and estimator must operate at high sample rates. Therefore, the estimator must be computationally 
efficient and capably process thousands of sensors at high sampling rates, assuming that all of the sensors along the 
fiber are utilized for feedback. Since more sensors leads to more efficient estimates, computational efficiency is a 
key feature of a robust modal filter. 

The estimator must not be influenced strongly by sensors located at leverage points. Leverage points are sensor 
locations which strongly influence feature estimates. These points are found on most structures and can be used to 
determine optimal sensor placements. 11 Leverage can be thought of as a moment arm. An outlier at a place having a 
large moment arm can drag OLS estimates significantly away from the majority of the data. 

The variation of the estimates between each discrete time step must be small, or the controller and system may 
become unstable. The selection of estimators is thus limited, as many estimators rely on random sub-sampling. 12 ' 13 
Approximate or stochastic algorithms may lead to inconsistent estimates in a non-convex optimization problem. 14 

Another significant challenge is that in the literature, there may not be a dedicated robust estimator to meet all of 
these requirements. A plethora of robust estimators exist for static analysis. A few are the Least Trimmed Squares 
(LTS), 12 Fast-S, 1 ' and Repeated Median (RM). 15 Robust estimators are often not required to be real-time estimators, 
as they find many uses in processing complex data analysis. Most robust estimators are computed over a period of a 
few seconds, minutes, or hours. Computation time often depends on the size of the data population. Most robust 
estimators assume initial unbiased population distributions, such as the normal distribution. 

In order to utilize robust regression methods, the estimator will have to perform near the controller sampling rate, 
which may be on the order of a few milliseconds. The majority of multivariate robust estimators are not real-time 
estimators. 14 

Some estimators smooth data over short time windows, where random outliers may occur, 16,17 but the bias 
induced by a fiber failure affects the modal estimates permanently after the failure. Thus, time window operations 
will not be sufficient, as the biased data become the majority of the data within the time window. 

The most promising estimators utilize two stages. The first stage produces a high breakdown estimate, and that 
estimate is refined in a Gaussian-efficient estimation stage. Modified Maximum Likelihood estimators 
(MM-estimators) initialize a Maximum Likelihood Type estimator (M-estimator) 18 with a feature estimate from a 
high breakdown estimator. 19 The M-estimate inherits the breakdown point from the original estimator and is as 
efficient as the final estimator. Relationships between M and S estimators have been formed 20 for MM estimators. 
Generalized M-estimation (GM) 21 has been developed, which down-weights leverage points. Down-weighting 
sensors due to high leverage leads to loss of efficiency, and is not pursued. Unfortunately, a major challenge is that 
two-stage estimators are impractical to compute because the initial estimate is still computationally intensive. 14 The 
next section presents an overview of the approach taken to overcome all of the aforementioned challenges. 

C. Overview of Approach 

In order to address the many challenges of an on-line robust distributed-sensor system, a new estimator called 
the concentrated modal estimator (CME) is derived. The CME is best described as a symbiotic estimator merging 
ideas from Tukey’s redescending M-estimator and concentration principles. Robust estimators based on 
concentration operators (CO) are used in the robust regression community 14,22 and have been shown to be 
computationally simple, consistent, and high breakdown. The M-type estimators are widely used and useful for 
many data distributions, because they can be tuned to be asymptotically normal and robust to most outliers. 

The redescending M-estimator is solvable through computationally efficient iterative recursive least squares 
(1RLS). The M-estimate uses weights which are tuned over IRLS iterations. The weights are a function of the scaled 
residuals, where proper scaling can lead to estimates which are 95% Gaussian-efficient. 18 But the M-estimator has a 
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0 breakdown point in the presence of biased sensor data at leverage points. Just as OLS, even one outlying sensor on 
a leverage point can drag the estimates away from a good solution. This fact is significant because the FOS system 
on a typical structure is full of leverage sensors. 

Concentration algorithms are only applicable to nominally multivariate normal data. The algorithms achieve 
high breakdown estimators by relying on iterative application to a smart selection of sensors with low breakdown 
estimators such as OLS. This method will not work in the presence of nominally asymmetric data. 

To shore up the limitations of both estimators, two approaches are taken. The M-estimator is initialized with a 
robust start (or initial feature estimate) 23 based on the previous time step estimates. This start is not susceptible to 
outliers. A robust trimming criterion in the concentration operator is also introduced. The criterion accounts for 
known asymmetric outliers and is derived from the Mahalanobis distance of the strain modal matrix data. Together 
with the M-estimates of location (mean) and dispersion (covariance), the trim criterion is applied to iteratively trim 
out bad sensor data. 

It is shown with Monte Carlo simulation (MCS) that the CME estimates converge for the aeroservoelastic trim 
case and for wing bending and torsional perturbations thereof. By comparing the CME to the state-of-the-art 
M-estimator with Huber or Tukey bisquare weightings, it is shown that the CME gives estimates with less bias for 
similar computational speed. Finally, the CME is demonstrated in dynamic ASC simulations on the X-56A model in 
the presence of 230 simultaneously failed sensors. For this scenario the robust modal filter and shape controller 
achieve ASC, but with expected tracking errors due to residual modes. The next section introduces the robust modal 
filter methodology which includes the derivation of the CME. 

II. Methodology 

The previous sections show the need for a modal filter to address the problem of outlying sensors on leverage 
points and that of computational efficiency. The following methodology is developed to directly handle this 
challenge. The CME is derived to robustly estimate the modal displacements of the X-56A aircraft (see Fig. 1) 
during a fiber optic sensor failure. The CME is a real-time concentration algorithm using robust starts, M-estimates 
in the concentration steps, and a fixed robustness trim criterion. The functional architecture of the CME is presented 
in Fig. 5. 
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outlying sensors 
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small outlying sensors 
Efficiency improvements 
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( ¥ c f ) 


Robust modal 
displacements 
at time x 


Figure 5. The functional architecture of the concentrated modal estimator. 

From Fig. 5, it is evident that there are two feedback loops. An inner loop represents the iteration of steps 
(called M-steps) of the M-estimator. The outer loop represents the trimming of gross outliers; each time through the 
loop is referred to as a concentration step. The sensor data flow into the system, and a first estimate of the weights 
and modal coordinates is initially computed. Outlier-sensitive information such as previous modal coordinates and 
weights are used to trim data in the concentration step. A new “good” sensor set is formed from the concentration 
step. This information is then passed to the M-estimator, where weights and modal displacement estimates are 
recomputed with iterative M-steps. The output of the M-estimator moves to the input of the concentration 
procedures; this process continues until convergence. This figure should be referred to as a guide for understanding 
the following sections. First, the strain-based M-estimator portion of the CME is derived. 
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A. Strain M-estimator Derivation 

The asymmetric nature of the distribution (see Fig. 4) demands a more robust estimator within the concentration 
procedure, such as the computationally efficient M-estimator. M-estimators are characteristically gradient descent 
algorithms. 2 They are computationally efficient, affine equivariant, robust to masking effects, and tend to 
outperform OLS when applied to many data sets. 25 Maronna’s Robust M-estimator 26 and a concentration 
algorithm 27 have performed similarly well for contaminated data sets used in principal component analysis." 5 
Merging the two concepts is attempted here. Indeed, Olive suggests that robust estimators can be used in place of the 
classical estimator for a concentration algorithm in some cases. 23 

The low theoretical maximum breakdown point of 1 / ( m + 1) of the M-estimator 29 is inconsequential for two 
reasons. First, this breakdown point is computed assuming that outliers can occur in all features of the data matrix. 
For the FOS system, outliers can only occur in one feature of the data matrix - the time-varying sensor measurement 
vector. Any outliers in the fixed portion of the data matrix will be accounted for with a trim criterion. The second 
reason is that a concentration operator does not require a high breakdown estimator within the concentration steps to 
lead to a high breakdown estimator. In fact, the concentration algorithms which employ OLS have been shown to be 
high breakdown for nominally multivariate normal distributions. 

The strain at measurement locations may be expanded as a summation of an infinite number of orthogonal strain 
mode shapes 30 as in Eq. (5). 


s k (x c ,y c ,z c ,t ) =2_ i q i (t)ip i (x c ,y c ,z c ) 


( 5 ) 


To reduce model complexity, only a subset m of mode shapes which dominate the response are generally included 
in the modal matrix. 3 It is assumed that the subset of modes captures the main dynamics and the sensors are subject 
to random errors. This error can be modeled as a normal distribution £ k £ N (ji n , a n ). At any discrete time step, 
t — t, the quasi-static approximate reading of any sensor can be given by Eq. (6): 


s k {x c ,y c ,z c 


.0 = ^ qi(T)ipi(x c ,y c ,z c ) + £ k 


( 6 ) 


where m is the number of mode shapes retained in the model. Consider the linear model for the k th sensor 
measurement to be described by Eq. (7): 


s k ( x c>yc> Z C> 


T) = ^ qi(j)lpi(.X c , 


yc.Zc ) + e k = 'V k (x c ,y c ,z c )q{ r) + e k 


( 7 ) 


where e k is a finite residual (that is, measurement error), v V k (x c , y c , z c ) £ M 1 xk is the k th row of the strain matrix, 
and q(j) £ M kxl is a vector of estimated modal displacements. From the sensors readings, the objective is to 
estimate q(r). This equation can be solved as a maximum likelihood estimation (MLE) problem (see Ref. 18) which 
is posed as minimization of an equally weighted summation of a function of the residuals as in Eq. (8): 

s s 

2^P(e fc ) =2^p(s k (x c ,y c ,z c ,T) -^(xcycZcjqfr)) (8) 

fc = i fc = i 


where S is the set of strain sensors and p(x) is an objective function with special properties. A reasonable p(x) must 
be even, zero when evaluated at zero, increasing for increasing arguments, and differentiable. Define the influence 
function cp(x ) = p’(x) as the differential of the objective function p(x). The influence function characterizes the 
proportional impact of the residuals on the estimate. The impact of an OLS residual on the estimate is directly 
proportional to the size of the residual, which is why OLS is not robust. To find q( r) the summation given in Eq. (8) 
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is differentiated by q(r) and is set equal to zero. By completing this operation, the equality shown in Eq. (9) is 
achieved. 


s 

^ <p{s k (x c ,y c ,z c , t) - ^ k (x c ,y c ,z c )q(r))H' fc (x c ,y c ,z c ) = 0 (9) 

k = 1 


Let w k (e k ) = C P(~ ek ^/e k for any <p(e k )- th en the weighted objective function can be rewritten as in Eq. (10): 
s 

^w fc (e k )(s fc (z c ,y c ,z c ,T)- l E k (z c ,y c ,z c )q(T)) l P k (z c ,y c ,z c ) = 0 (10) 

k = 1 


which results in the weighted least-squares problem. 31 For all sensors, Eq. (10) forms a system of equations, which 
when solved give an efficient estimate of q(r) under normal conditions. The weights w k (e k ) are affine equivariant 
and modeled as functions of the residuals, e k , and the residuals are functions of the weights. Therefore, recursion 
(that is, IRLS) is required. This operation proceeds by solving for an initial least-squares estimate q(r) and 
computing the residuals and weights. Using the weighted observations, a new feature estimate q(r) is computed, 
and the residuals and weights are recalculated. The features or modal displacements q(r) of the hyperplane 
approximately satisfying for all sensors, Eq. (10), usually appear within a few iterations. 

B. Weight Computations of Strain M-estimator 

The solution of Eq. (10) must be computed after each concentration step, c, for the proposed concentrated 
estimator (see M-steps in Fig. 5). To improve the convergence to the unbiased solution of q(r), sensors which are 
the most outlying are completely removed. For the new group of sensors, M-estimation is used to find improved 
feature estimates. Selection of the influence function is critical to the performance of the M-estimator. 

Two commonly used influence functions in M-estimation are the Huber function 18 and Tukey’s bisquare 
function.'" While robust and efficient in many cases, Huber’s influence function increases without bound as the 
residual departs from 0. Therefore, gross outliers still impact the feature estimates and in typical cases lead to 
efficiency losses of 10-20%. 33 

Tukey’s bisquare function belongs to a class of redescending functions 34 which account for gross outliers by 
gradually reducing the influence of the large residuals. Redescending M-estimators use (p{x) influence functions 
which are non-decreasing near the origin, but decrease to 0 far from the origin at some minimum rejection point. 

For this reason, Tukey’s bisquare function is chosen to compute the weights with the residuals of the data. The 

bisquare weighting function w = ( P^ X ^/. % is defined for the k th sensor as in Eq. (11), 



, 0 otherwise 


where is the median absolute deviation (MAD), h 0 is a tuning constant, c is a concentration step, and b is an 

IRLS step of the M-estimator referred to as an M-step. To achieve the maximum 95% asymptotic efficiency 
assuming residuals have a Gaussian distribution, it has been shown that a tuning constant of h 0 — 4.685 is 
required. 35 The MAD for the k th observation is calculated as in Eq. (12): 

(6 , c) _MED{\e^-MED{e^)\)l 

a k ~ / .6745 (12) 
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where the constant scaling 0.6745 is required to achieve a 37% Gaussian-efficient consistent estimator of the 
standard absolute deviation. ’' 1 While relatively low-efficiency, the purpose of using MAD instead of using the true 
scale is to resist outliers. This resistance it achieves remarkably well, because the median is high breakdown. 
The MAD is developed for symmetric distributions and does not address distribution skewness, which may be of 
concern since the explanatory data (strain mode matrix) is multivariate-skewed. Improvements of the MAD 
approximation for asymmetric long-tailed distributions are available if necessary (see two alternatives in Ref. 36). 
Given the weights, w ( k b,c \ the linear system of equations is solved for q^ b,c \ t), given sensors in subset Sg as in 
Eq. (13). 


) ( _ ^))( s fc^c-yc.z c ,T) -T' k (x c ,y c ,z c )q (b ' e) (r))T' k (x c ,y c ,z c ) = 0 ( 13 ) 

/c = i / 

The weighted least-squares problem for the c th concentration step is solved in the same way as in Eq. (10). 
Equations (11)-(13) are the primary feature estimator equations used within a concentration step of the concentration 
operator. These equations are iterated within any concentration step for a specified number of M- steps, bf resulting 

in the c th feature or modal displacement estimate, q^ b f' c \x). 

C. Concentration Operations 

The output of the M-estimator is prepared for concentration in Fig. 5. The purpose of concentration is to 
iteratively remove poor observations (sensor measurements) and use the sensors that are closest to the unbiased 
estimate of the centroid of the data distribution. Utilizing the sensors nearest to this centroid is assumed to give the 
best feature estimates. A best estimate of the centroid is the multivariate location, , and dispersion, V , of the data. 
Redescending M-estimators have been proposed as robust estimators of multivariate location and dispersion for 
theoretical asymmetric distributions. ,7 

Sensors furthest from this centroid are proposed to be downweighted in Eq. (13); however, downweighting 
sensors puts initial trust in possibly gross outliers. Therefore, the most offending observations must be completely 
removed from consideration. 27 Although the redescending M-estimator does in fact equate weights to 0 for gross 
outliers, it puts some initial trust in gross leverage outliers in the first M-step. It is shown later that converged feature 
estimates from a redescending M-estimator may remain biased in some cases due to leverage outliers. The weighted 
sensor removal methodology to improve gross outlier rejection is developed here. Let the k th sensor data vector be 
defined as in Eq. (14). 


X k = [v k (x c ,y c ,z c ) s k (x c ,y c ,z c , t)] (14) 

Defining the data row vector in this way ensures that time-varying outliers shall not occur in m features of the data 
matrix, because the strain modal matrix is a fixed (not dependent on t) and known data set, assuming no adaptation 
is present. Defining the data vectors in this way dramatically increases the breakdown point of an M-estimator. From 
any sample sensor set Sg 2 S, a location vector (see Eq. (15)) 






(15) 


and a dispersion matrix (see Eq. (16)) 
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( 16 ) 


v <^) = ^ w M x r Xk 


are estimated in the c th concentration step. The weightings are the result of br iterative M-steps over the subset of 
sensors Sg. Weighted location and dispersion matrices have led to robust equivariant estimators with a high 
breakdown point for any dimension, 38 such as the Stahel 39 and Donoho 40 estimator. It was shown that if the weights 
are affine equivariant, the estimates of location and dispersion are also affine equivariant. It was also shown that if 
the true mean and dispersion of the model has an asymptotic breakdown of 0.5, then the asymptotic breakdown 
point of the location and dispersion estimates also have an asymptotic breakdown of 0.5. 

For the estimated location and variance, the squared Mahalanobis distance ( D 2 ) (see Ref. 41) is computed for 
every sensor data point k by Eq. (17). 

D 2 (x k ) = (: x k - T (fe /' c) ) 1 ( x k - T&frf (17) 

This multivariate distance differs from the Euclidean distance only in that it accounts for correlations between data 
points and is scale-invariant. If the population has a multivariate normal distribution, the D 2 is asymptotically 
approximated by a chi-square distribution. 8 With this knowledge, statistical cutoff points from the inverse 
cumulative distribution can be determined. The strain data matrix has, however, an unknown highly skewed 
distribution, thus this data removal technique will not succeed. 42 Theory-based concentration algorithms which trim 
the percentage of observations having the highest D 2 are strictly invalidated. Leverage points naturally have very 
large D 2 , therefore, trimming good leverage points drastically biases the feature estimates. 

The amplitude of D 2 remains useful for finding outliers if the multivariate normal assumption is violated; 
however, asymptotic theoretical cutoffs must not be relied upon. Without knowledge of the underlying theoretical 
distribution, an approximation is required to find the cutoff value of D 2 . The initial distribution of D 2 may be 
computed from the fixed modal matrix and time-varying set of strain data with Gaussian noise. The maximum of 
the computed D 2 may be used as an upper bound for removing gross outliers. This method is very similar to the 
empirical cutoff approach for a fixed data set described in Ref. 43; that approach was improved with the adaptive 
approach taken in Ref. 42. 

A shortcoming of these two methods is that small outliers may be missed if sensors are removed based on a 
maximum threshold of D 2 or some derivative method, because the initial distribution mean and covariance may be 
biased. Iterative concentration steps are proposed herein to address this problem. During each concentration step, 
gross outliers are removed and the location and dispersion are re-estimated. The sample location and dispersion 
more closely resemble the population location and dispersion, thus, the small outliers become more pronounced. As 
the D 2 (x k ) increases, the sensor can be identified as an outlier and removed. Outliers missed by this trim procedure 
are more likely to be down-weighted in the M-estimate (see Eq. (13)). 

The proposed method for finding the upper bound D 2 b is time-consuming to implement, requiring thousands 
of simulations because the strain is time-varying. Since most of the data are described by the constant-strain 
data matrix, an approximation can be used for the upper bound. It can be assumed that the distribution of 
D 2 (x k ),k = 1 ...5 is equal to or greater than the distribution of D 2 { l V k (x c ,y c ,z c ) s ),k = 1 ...S if the sensor data 
have a Gaussian error distribution. With this assumption, the impact of an additional feature may be assumed to 
change the distribution of D 2 by the additional degree of freedom impact in a chi-square distribution. Recall that D 2 
is given in units of variance, which implies that the variance will increase with the additional degree of freedom. 
Therefore, it can be assumed that D 2 ( (x Cl y c , z c ) ) + $(n s ) > D 2 (x k ), k = 1 ...S. Assuming the adjustment of 
$(n s ) is due to the noise of the strain data, the scalar upper bound is defined as shown in Eq. (18): 

D ub = Pc max D 2 ( V k (x c , y c ,z c )) ( 1 8 ) 

where P c is a tuning constant chosen to be slightly greater than 1. The tuning constant accounts for i9(n s ). By 
removing a portion k £ S sensors with D 2 < D 2 b , a new candidate group of sensors S£ +1 is found for the next 
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concentration step, and consecutive M-steps. The bound proposed in Eq. (18) is both theoretical and empirical and is 
the meat of the concentration procedure (see Fig. 5). Simulation studies given later verify this approximation of the 
upper bound, D/ h to be good for the strain mode matrix and strain data. The next section discusses another method 
of improving the robustness of the concentration estimator. 

D. Robust Starts and Online Operations of the Concentration Modal Estimator 

Robustness for multi-stage estimators tends to come from good starts (initial feature estimates). A feature 
estimate from a high breakdown estimator is used to start the M-estimator for MM-estimates. 14 The robustness is 
inherited by the more efficient M-estimator; however, this process can be time consuming because most high 
breakdown estimators are computationally inefficient. This attribute presents a problem for a distributed-sensor 
system, which requires a high breakdown estimate but must also be computationally efficient. 

Other concentration operators use starts from estimates from all of the data or the data closest in distance to the 
coordinate-wise median of the data. The median ball algorithm 23 uses feature estimates from sensors closest to the 
median as a robust start. This is a good start if the data can be assumed to be multivariate-normal, and works 
reasonably well for skewed distributions. 

The first estimate of the system when t is 0, (that is, when the sensor system is first operational), is calculated 
with a non-robust least-squares estimate. The first estimate is assumed to come from a working sensor system, thus 
it is a robust estimate. The initial robust feature estimate q fo ' oi (0) is found by solving the least-squares problem 
shown in Eq. (19): 


^ (s fc (x c ,y c ,z c ,T) -H , k (x c ,y c ,z c )qC°'°)(0)) l F k (x c ,y c ,z c ) = 0 

k= 1 


(19) 


where Sg is the set all of the available working sensors. 

During operation, a robust start is paramount. A significant advantage of a time -based sensor system is that 
previous close estimates are available. The most robust start will be the estimate from the previous time step, 
because the strain change is expected to be small between discrete time steps. Thus, the robust starts between 
discrete time steps are implemented as shown in Eq. (20): 

q (0 ' 0) (r) = q ib f' c f\T- 1) (20) 

where bf is the total number of M-steps chosen, and cy is the total number of concentration steps. 

The importance of starts carries over into the concentration steps themselves. In order for the steps to be high 
breakdown, each concentration step requires a robust start. The initial start, q f0,Hi (0), being robust, the final 
estimates at the end of each of the concentration steps: ...q^ b f' c f~ 1 \T') are robust under the 

assumption of robust inheritance. 14 The estimates of corresponding concentration steps, then, are robust starts for 
the respective next concentration steps: q^ 0,2 ^(r), q(°' 3 )(r), ... q(°<<y)( T ) Therefore, the inheritance assumption 
shown in Eq. (21) is used to generate robust starts between concentration steps (see outer-loop feedback in Fig. 5): 

q(°' c+1 >( T ) = q (h f' c \r) (21) 

The full steps of the CME for any discrete time step are summarized in Algorithm 1 (also see Fig. 5), assuming that 
an initial feature estimate has already been calculated with Eq. (19) at time 0. 

Algorithm 1: {s k (x c ,y c ,z c ,T),q^ b f' c ^(j — 1)} -> q^ b f ,c f\ t) 

1. Ifc = 0, compute q (0,0) (r) using Eq. (20); otherwise compute q (0,e) (r) using Eq. (21). 

2. For b = 0: bf, iteratively compute weights, w^ bf ' C \ using Eqs. (1 1)-(13) to obtain w^ f ' C \ 
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3 . 


Compute location T (see Eq. (15)) and dispersion V (see Eq.(16)) with w^ b/ ' C \ 

4. Compute D 2 (x k ),k = 1 ...S, using Eq. (17) with T and S. 

5. Generate a new sensor set 5^ +1 by trimming sensors below cutoff D 2 h in Eq.(18). 

6. If c < Cf, go to step 1; otherwise output q( 6 /' c /)(r). 


For each time step, the M-step iteration count bf may be initialized to be large, so that a robust redescending 
M-estimate initializes the CME. This improves the algorithm’s stability during the concentration steps. Afterward, 
single M- steps where h f is equal to 1 may be utilized. This method has the effect of improving computational 
efficiency. 


E. Concentrated Modal Estimator Summary 

The CME finds great application for sensor systems with very large numbers of data points, such as will be 
available with the FOS, because some of the sensors are just not as important as others and there are hundreds from 
which to choose. If some sensors are biased, others can be used in place of those biased to form modal estimates. 

The CME is noticeably similar to previously derived estimators. The CME uses concentration steps as proposed 
for the DGK (Devlin, Gnanadesikan, and Kettenring) estimator (see Ref. 14) and median ball algorithm proposed by 
Olive. 23 Rather than removing a percentage of data at every concentration step, however, a datum is only trimmed if 
its D 2 exceeds D 2 b (see Eq. (18) and step 5 of Algorithm 1). The estimator thus follows the Hippocratic Oath, which 
may be paraphrased as “do no harm.” 

Another difference includes robust start inheritance used between concentration steps (see Eq. (21)) and between 
time steps (see Eq. (20)). The median ball algorithm uses two starts, including the median start and the classical 
start, because access to close estimates of population parameters is not available. A previous close sample estimate 
will likely outperform a geometrically robust start, especially if the data are heavily skewed. 

The CME is a deterministic algorithm, and requires no random subsampling. Most robust estimators rely on 
random subsampling; an example is the popular LTS estimator." 9 It has been shown, however, that estimators with 
random seeds are not consistent. 14 Instability may result if large (incorrect) changes in modal displacement estimates 
occur between time steps. The deterministic approach of the CME leads to stable estimates which do not vary by 
re-running the algorithm. 

It is difficult to see how the deterministic concentration procedure or the start can negatively affect the 
high-breakdown nature of the redescending M-estimate. With robust starts and high-breakdown implications over 
time and over concentration steps, robustness will likely be achieved by the CME. In fact, the breakdown point can 
be higher than 0.5 due to the robust start utilized in Eq. (20). Simulation studies presented later justify the CME as a 
robust estimator for several worst-case asymmetric data distributions. 


III. Simulation 

The CME is demonstrated in static and dynamic simulation studies. First, the sensor failure simulation is 
developed. An appropriate worse-case scenario failure point is determined. For analyses, a failure in a fiber is 
induced in a critical location. The CME is applied to scenarios in which the wing is in aeroservoelastic trim and 
perturbed from trim. An aeroservoelastic trim includes trim modal displacements. 

Monte Carlo simulation is used to gather error distribution estimates for the modal displacement approximation. 
The modal estimate errors are compared to the state-of-the-art M-estimator feature estimates. A computational time 
study is given to show the CME has the potential to be a real-time estimator. Finally, a dynamic simulation verifies 
that the ASC system for the X-56A model does not go unstable. This includes a comparison of the CME to the 
state-of-the-art robust estimator. 


A. Fiber Optic Sensor Failure Simulation 

The controller requires accurate modal estimation to track the displacements at the locations given in Fig. 2. 
Modeling FOS failures is required to test the robustness of the modal filter and FOS-based control system. 
Researchers at LaRC investigated the nature of spurious strain data after a break in the FOS fiber occurred. 4 From 
visual inspection of the data it appeared that high bias in the strain occurs just upstream of the break in the fiber. 
Downstream of the break, the strain measurements appear biased to have a mean of 0 and a low standard deviation. 
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These characteristics are captured here for a SFOS failure; however, this particular failure model may not be the 
general case. This demonstration should, however, lead to systems which model small, medium, and gross outliers. 

The sensor locations P*f(x,y,z) upstream (closer to the wing root) from the fault location Pf (x, y, z) are found, 
within a radius, r n f. The relative bias shape on the k th sensor upstream of the fault is modeled by a normal 
distribution as shown in Eq. (22). 



r n/ V27T 


exp 


1 { \\Pf(Xc,yc,Z C ) - Pnf( x C,yc,Z C )\\ ' 


2 V 


n f 


( 22 ) 


The sensors nearest the fault are modeled to have the most bias; those farthest from the fault are modeled to have the 
least bias. The bias is added to the sensor measurements with the rule shown in Eq.(23): 


S^ f (x c , y c , z c , t) = (x c , y c , Z C , t) + 


- c n f ( 


B 


nf 


max B 

keS nf 


nf 


A 


(23) 


where A is the maximum desired strain variation on sensors upstream of the fault in the fault radius. The sensors 
downstream of the fault (outboard near the wing tip in this case) also experience a bias; however, rather than a bias 
added to the existing measurement, the bias is modeled to take over the sensor measurement completely. This 
condition is modeled by replacing the sensor measurement with a sample from a normal distribution with a mean of 
0 and a standard deviation of half the magnitude of A, as shown in Eq. (24). 

Sk f ix c ,y c ,Zc,r) = N( 0,^) (24) 

The amplitude is divided by 2 to make the error variation smaller farther from the fault. The mean was selected 
to be 0, but this may vary depending on the way the fiber is failed. Certainly this is not a perfect model of a fiber 
optic sensor fault, but the bias added to the sensors using Eqs. (22)-(24) is appropriate for demonstrating outlier 
rejection. 

For the failure model given in Eqs. (22)-(24), three structural strain scenarios are analyzed. The first structural 
strain scenario is for aeroservoelastic trim strain at the design speed. This is a strain scenario in which the aircraft 
wing will spend the most time. The second structural strain scenario is for a large wing tip leading-edge -down 
torsional displacement from aeroservoelastic trim. The third structural strain scenario is for a large-amplitude 
bending displacement from aeroservoelastic trim. 

It is expected that large displacements from aeroserovelastic trim may result from maneuvers, shape control, or 
large disturbances. To simulate the expected failure bias during a break, the failure bias amplitude, A , is arbitrarily 
set to 30 times the standard deviation of the SFOS noise (see Eqs. (22)-(24)). The SFOS normal error was assumed 
to be 3 microstrains (jis) because the FOS is expected to have a high signal-to-noise ratio. The radius, r n f , which is 
used to find biased sensors upstream of the fault, is set to 3 inches. 

The radius selection is somewhat insignificant, as the break in the SFOS is assumed to occur near the trailing 
edge of the right wing near the wing root. This location has the highest leverage points or Mahalanobis distance 
(see Fig. 4). The nominal sensor measurements for all three scenarios superimposed with suitable sensor bias for a 
fiber break at the wing root is presented alongside the other five fibers in Fig. 6. 
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Figure 6, Simulated fiber optic sensors strain with fault + noise: a) trim strain; b) trim + torsional strain; and c) trim + 
bending strain. 

The biased strain in Fig. 6(a) represents small outliers. The biased strain in Fig. 6(b) caused medium outliers. 
The biased strain in Fig. 6(c) represents a case with gross outliers which occur both at non-leverage and leverage 
points (see Fig. 4). The strongly biased strain measurement data (see Fig. 6) present a challenge to the modal 
filtering and control system. 

B. Concentrated Modal Estimator Simulation and Comparison 

For each structural strain scenario. Algorithm 1 is computed for four concentration steps. The CME requires a 
robust start from a previous time step; in this case the robust start was not available. The robust start was therefore 
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modeled by the true modal displacements offset by 10% of a multiplicative normal error. The relatively large offset 
simulates the modal displacement variation between discrete time steps r. Recall that modal displacement estimates 
current discrete time steps are used as robust starts for future time steps in the CME. The number of M-steps in c 0 is 
initially set to 10 to achieve a converged Tukey bisquare M-estimate and then is set to 1 for all remaining 
concentration steps c 0 , c 1( ... , cy to improve computational efficiency. The tuning constant P c for the D/ h required for 
each concentration step is set to 1.1. The D^ b works out to be 68 using Eq. (18). 

The CME is compared to M-estimates with Huber and Tukey bisquare weightings. Huber’s function is utilized 
because it down-weights but does not completely remove the presence of gross outliers. Its performance is 
comparable to that of OLS used by Kang et al; 44 however, it will be much more robust to outliers. The M-estimators 
are given the same robust start as the CME: the true modal solution offset by 10% multiplicative normal noise. 
Recall that the additional noise simulates the difference in modal estimates between time steps. 

Since control systems require high sampling rates, the CME must have low computational complexity. The 
computational processing time used for all estimators is recorded with the MATLAB profiler, which estimates the 
total CPU time required by processors to run functions and sub -functions. For each scenario a 2.6-GHz processor is 
used to compute CPU time. Since the noise and fault conditions are characterized by normal distributions, an MCS 
is run. The MCS is generated from 300 random seeds. 

Results are presented for percent relative error and deviation for modal bending and torsion displacement 
estimates. The simulation modal displacement is considered the true model of modal displacement in the system. 
The results of the MCS simulations for the aeroservoelastic trim case and the structurally perturbed cases are 
presented in Fig. 7. 
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Figure 7. Modal estimates during fault: Relative error a) trim strain, b) trim + torsional strain, and c) trim + bending 
strain; and CPU time d) trim strain, e) trim + torsional Strain, and f) trim + bending strain. 

The significance of Fig. 7 is primarily in the relative error comparisons. In the first structural strain scenario in 
Fig. 7(a), the relative error distribution of the SW1B modal displacement estimated with Fluber weights is 
symmetrical and centered at -7%. The first standard deviation moves the overall maximum error to -15%. The 
SW1T modal displacement relative error distribution is skewed negatively and centered at -3%. The maximum 
deviation of the error moves the error to -22%. 

The Tukey estimates fared better, but only slightly. The SW1B modal displacement error distribution estimated 
with Tukey functions is symmetrical and centered at 0%. The error deviation is up to 10%. The SW1T modal 
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displacement error distribution is symmetrical and centered at 2%. The maximum deviation of the estimate goes up 
to 18%. Reduced mean errors are expected for Tukey function estimates due to the reduction of the influence of 
gross outliers with bounded influence functions. The error bars were nearly the same size for both estimators. 

The CME estimates the SW1B modal displacement with an error distribution for both SW1B and SW1T modal 
displacements symmetrically centered at 0.5% in Fig. 7(a). The deviation of the error for the SW1B modal 
displacement was at a maximum of 5%. The deviation of the error for the SW1T modal displacement achieved a 
maximum of 10%. When compared to state-of-the-art estimates, CME outperforms them with respect to relative 
error for the aeroservoelastic trim case. Figure 7(d) indicates that the CME is computationally comparable to the 
state-of-the-art estimators. The means of the CPU time for the CME was at 25 ms. The CPU time varied 18 ms from 
the mean. 

Fig. 7(b) shows the relative error comparisons for the second scenario, in which the wing is elastically twisted 
leading-edge-down by 3 deg. With higher displacements from trim, the estimators are expected to perform worse, 
due to the growth of outliers, and in fact this is shown to be the case. The Huber SW1B modal displacement error 
distribution is skewed positively and centered at 7%. The maximum deviation of the error moves the relative 
error up to 14%. The SW1T modal displacement error distribution is skewed negatively and centered at -28%. The 
error variation takes the maximum error to -47%. 

Tukey’s estimate is better than Huber’s but worse than for the aeroservoelastic trim scenario. The SW1B modal 
displacement error distribution is skewed negatively and centered at 3%. The maximum relative error is down 
to -10%. The SW1T modal displacement error distribution is symmetrical and centered at -7%. The error variation 
takes the error distribution to -20%. 

The CME estimates for the torsional scenario are comparable to the aeroservoelastic trim case. The means of 
both modal estimates are symmetrical and centered near 0%. The SW1B modal displacement estimate varies up to 
4% in either direction. The SW1T modal displacement distribution varies up to 8%. The CME outperforms both the 
Tukey and Huber estimates. The CPU time for the three estimators shown in Fig. 7(e) is nearly the same as for 
the aeroservoelastic trim case, however, the CME CPU time distribution increased to 31 ms with a 17-ms variation. 

In the final scenario the biggest improvement is seen when using the CME compared to the Huber and Tukey 
estimates. Huber’s estimate is strongly biased. The SW1B modal displacement error distribution is nearly a point 
and centered at 7%. The SW1T modal displacement error distribution is symmetrical and centered at 145%. The 
error varies up to 190%. The torsional modal displacement estimate is extremely poor. This holds true for the Tukey 
estimate as well, the SW1T modal displacement distribution of which is symmetrical and centered at 20%. The error 
variation of the estimate is up to 48%. 

The CME estimate shows almost no error bias in the SW1B modal displacement. The SW1T modal 
displacement error distribution is higher than from previous scenarios, however, the mean is near 0 again. The 
variation is up to 20%. The clear advantage seen in the third scenario comes from the handling of gross outliers at 
leverage points through the concentration procedure. Neither the redescending M-estimator based on Tukey’s 
bisquare function nor the M-estimator with Huber weights considered significant removal of these leverage outliers. 

C. Analysis of Concentration Steps 

The previous results are telling of how the CME will outperform the state-of-the-art estimators for the 
asymmetrical multivariate estimation problem. The CME process of concentration is not completely intuitive 
without analysis of the squared Mahalanobis distance D 2 at each concentration step. For the aeroservoelastic trim 
strain scenario, the initial distribution of D 2 is given, along with the measured D 2 and weighted D 2 for four 
concentration steps. The initial distribution is based on D 2 ('¥ k (x c ,y c ,z c ')'), k — 1 ... S, where the squared 
Mahalanobis distance is computed only for the fixed-strain mode matrix. The measured D 2 includes the strain mode 
matrix and measured strain in the computation of the squared Mahalanobis distance. The weighted D 2 is computed 

by multiplying the measured D by the final weights w fc from the CME for each sensor. For the aeroservoelastic 
trim strain scenario, the distribution of the D 2 for all SFOS is given for successive concentration steps in Fig. 8. 
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Figure 8. Mahalanobis distances at the end of each concentration step in the aeroservoelastic trim case: a) c 0 ; b) c x ; c) c 2 ; 
d) c 3 ; and e) c 4 . 

Figure 8 gives several indicators that the CME is operating as predicted during its derivation. The first indication 
is that the measured and weighted D 2 tends to decrease through further concentration. At the beginning of the 
concentration procedure (see Fig. 8(a)), the measured D 2 is very large - up to 6,600 - and largest where the sensors 
have initially failed. The second concentration step, in Fig. 8(b), shows that the magnitude of the weighted and 
measured D 2 has reduced to a maximum of 250. In the final step (see Fig. 8(e)), the D 2 (x k ) of each sensor is below 
the D 2 b of 68. 

All of the sensors cannot be detected and trimmed in the first step because the mean and co -variance estimates 
are still biased. As the more biased leverage sensors are removed, the estimates move closer to the true population 
mean and covariance. As the true population mean and covariance are approached, the sensors with smaller bias 
begin to look more like outliers and cross the D 2 b threshold. These sensors are detected and removed, thus further 
improving the estimate of the mean and covariance of the distribution. This process is iterative and converging. 
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Notice from Fig. 8 that not all of the sensors can be removed with trimming, as outliers at off-leverage points 
are likely to reside below the D 2 b threshold. The effects of these outliers are down-weighted by the M-step 
re-weighting procedure. Since the weighted D 2 is below that of the good leverage points, the effects of these outliers 
have a minimal impact on the estimate. Therefore, the optimal feature estimates are pulled toward the true global 
optimum. 

Some computational observations of theoretical predictions can be made. Note that the measured D 2 is 
lower-bounded by the initial D z , supporting the fact that the addition of another feature and sensor noise to the 
initial D 2 increases the maximum D 2 . The utilization of Eq. (18) to approximate D 2 b is thus justified; it is best 
depicted in the last concentration step (see Fig 8(e), where the resolution is more pronounced. Another observation 
can be made about the effect of the weights on the noise: It is clear that the CME has a side effect of 
down-weighting noisy sensors; the weighted D 2 appears smoother than the measured D 2 . For those sensors which 
were particularly impacted by noise, the weighted D 2 was even lower than the initial D 2 . Thus, sensors with more 
noisy measurements than others can be identified and down-weighted within a single time step. 


D. Dynamic Simulation - Automatic Sensor Failure Rejection 

The previous static analyses show that the CME can perform adequately in the presence of unbiased and biased 
sensor data. But performance in a control system is a critical requirement of the CME, thus, the CME is tested in a 
dynamic simulation to verify that the interaction of the estimators with the control system will not lead to instability. 
For this verification test, the virtual deformation simulation control architecture shown in Fig. 9 is used (also see 
Ref. 2). 


Modal transformation 

► 



Figure 9. Virtual deformation control architecture for the X-56A model. 

Virtual deformation control is a concept which has been proposed to actively control the shape of the aircraft, 2 ’ 
however, to enable this concept in the manner previously proposed, distributed sensing is required. The distributed 
sensing and control architecture represents the control system for the simulated X-56A model, where the inputs are 
assumed to originate from a guidance computer. The commands are split into deformation and airframe type and 
the entire simulation and controller is run at 100 Flz. This sampling rate is faster than the predicted performance of 
the CME, but the algorithm has not yet been optimized computationally and placed into hardware. 

The simulated virtual deformation control system is thoroughly described in Ref. 2. For the present simulation, 
rigid body pitch angle, 6, and bank angle, 0, are tracked in the flight controller. Yaw axis commands are not given. 
Commands of 0 deg are given to both pitch and bank angles. Points at the wing tips (see Fig. 2) are given an equal 
1.2% of wing span vertical displacement commands. The points are tracked by commanding the first bending and 
torsion modal displacements using the transformation given in Ref. 2, explaining the “virtual” in “virtual 
deformation control.” 
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In previous work, 2 the simulation incorporated airframe noise, n AF , into the rigid body sensors. Only SFOS 
noise, n s , is modeled here, so that the effect of the fault is isolated. For the current simulation, the SFOS failure bias 
rif is added to faulty sensors using the same failure shown in Fig. 6(b). At any time after 10 s the sensor bias, ry, 
impacts the sensor system. The CME is allowed four concentration steps. As before, the CME is allowed 10 M - steps 
in the initial concentration step. A single M-step is utilized in the last 2-4 concentration steps. The D 2 h is again 
calculated to be 68, with P c set to 1.1 using Eq. (18). 

For comparison, simulation results for state-of-the-art Af-estimator with Tukey bisquare weights was utilized in 
lieu of the true state-of-the-art OLS estimator from Kang et al. 44 Clearly, an OLS estimator is an unfair comparison 
in the presence of such large sensor bias (see Fig. 6). During each simulation, both the CME and Tukey’s estimator 
use the robust starts in Eqs. (20) and (21). This is done to ensure a fair comparison and to demonstrate the 
importance of concentration. The Tukey estimator is also allowed to iterate to convergence. 

The modal displacement estimates from each estimate are passed to a /r-optimal controller developed in previous 
work. 2 The controller achieves robust stability and performance for modeled feature and speed variations there, 
however, it has some nominal overshoot performance problems, which may be corrected with improved weightings. 

It is not expected that nominal stability or performance problems from the controller will create other problems, 
thus, if instability occurs during the fault, it would not be expected to be the result of an improperly designed control 
system. Good or bad performance is due only to the estimator performance. The comparative results of the dynamic 
simulation studies are presented in Fig. 10. 
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Figure 10. Dynamic simulation comparing estimators during SFOS failure: M-estimate a) deformation tracking and b) 
airframe state tracking; and CME c) deformation tracking and d) airframe state tracking. 

Differences are noted from the side-by-side comparison of the M-estimator (see Figs. 10(a) and 10(b)) and CME 
(see Figs. 10(c) and 10(d)) performance. After 10 s, the control system with the state-of-the-art M - type estimator 
experiences strong divergent oscillations (see Figs. 10(a) and 10(b)), never returning to normal, and the system goes 
unstable. It is evident that the bias modeled by Eqs. (22)-(24) appears to lead to either control-induced instability or 
flutter amplification. This condition is worsened by the fact that the X-56A plant model is open-loop unstable. This 
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observation exposes the danger that may result from using a failed FOS system with an estimator which is not robust 
to leverage outliers. 

The CME time histories (see Figs. 10(c) and 10(d) show no signs of growing oscillations after the fault. When the 
structure is perturbed, the distribution of the noise does not appear to change and the estimates remain unbiased. The 
same is true when the structural command is changed, showing that the CME is consistently rejecting the outliers for 
different structural conditions. Note that “different structural conditions” also means “different outlier 
characteristics.” 

The dynamic performance of the CME is adequate when considering that 230 sensors have become strongly 
biased (as in Fig. 6(c)). The dynamic simulation demonstrates that the robust start between discrete time steps 
(see Eq.(20)) is justified in Algorithm 1 - that is, that the previous modal displacement estimate can be satisfactorily 
used as a robust start for the CME. 

The results presented here clearly show that the CME can help to create a robust distributed-sensor-based control 
system. Further simulation work may be necessary to verify that the CME is robust to all types of FOS failure modes 
in experimental studies. 


IV. Conclusions and Future Work 

The concentrated modal estimator (CME) was introduced as a candidate estimator to mitigate the severity of a 
fiber optic sensor failure. The algorithm works primarily under the assumption that there are many available sensors. 
The concept is that many sensors will work properly, and these can be used to determine which sensors are bad. 

The CME provides unbiased modal displacement estimates to the control system under simulated fiber optic 
sensor failure conditions. This estimation was achieved through the concentration procedure, which utilized an 
approximated squared Mahalanobis distance trim criterion. The CME was found to outperform state-of-the-art 
M-estimators, using the same robust starts. The CME was also shown to be computationally efficient relative to 
state-of-the-art M-estimators. 

The CME supports the safety-critical aspect of employing fiber optic sensors in distributed-sensing-based control 
systems. There are other applications of the CME, such as health monitoring and load safety, which are planned to 
be presented at a later date. Future work is also planned to extend to actively optimizing the shape of the aircraft. 
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